home *** CD-ROM | disk | FTP | other *** search
/ Sprite 1984 - 1993 / Sprite 1984 - 1993.iso / src / lib / m / cbrt.c < prev    next >
Encoding:
C/C++ Source or Header  |  1988-07-11  |  2.9 KB  |  102 lines
  1. /*
  2.  * Copyright (c) 1985 Regents of the University of California.
  3.  * All rights reserved.
  4.  *
  5.  * Redistribution and use in source and binary forms are permitted
  6.  * provided that this notice is preserved and that due credit is given
  7.  * to the University of California at Berkeley. The name of the University
  8.  * may not be used to endorse or promote products derived from this
  9.  * software without specific prior written permission. This software
  10.  * is provided ``as is'' without express or implied warranty.
  11.  *
  12.  * All recipients should regard themselves as participants in an ongoing
  13.  * research project and hence should feel obligated to report their
  14.  * experiences (good or bad) with these elementary function codes, using
  15.  * the sendbug(8) program, to the authors.
  16.  */
  17.  
  18. #ifndef lint
  19. static char sccsid[] = "@(#)cbrt.c    5.3 (Berkeley) 4/29/88";
  20. #endif /* not lint */
  21.  
  22. /* kahan's cube root (53 bits IEEE double precision)
  23.  * for IEEE machines only
  24.  * coded in C by K.C. Ng, 4/30/85
  25.  *
  26.  * Accuracy:
  27.  *    better than 0.667 ulps according to an error analysis. Maximum
  28.  * error observed was 0.666 ulps in an 1,000,000 random arguments test.
  29.  *
  30.  * Warning: this code is semi machine dependent; the ordering of words in
  31.  * a floating point number must be known in advance. I assume that the
  32.  * long interger at the address of a floating point number will be the
  33.  * leading 32 bits of that floating point number (i.e., sign, exponent,
  34.  * and the 20 most significant bits).
  35.  * On a National machine, it has different ordering; therefore, this code 
  36.  * must be compiled with flag -DNATIONAL. 
  37.  */
  38. #if !defined(vax)&&!defined(tahoe)
  39.  
  40. static unsigned long B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */
  41.                  B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */
  42. static double
  43.         C= 19./35.,
  44.         D= -864./1225.,
  45.         E= 99./70.,
  46.         F= 45./28.,
  47.         G= 5./14.;
  48.  
  49. double cbrt(x) 
  50. double x;
  51. {
  52.     double r,s,t=0.0,w;
  53.     unsigned long *px = (unsigned long *) &x,
  54.                   *pt = (unsigned long *) &t,
  55.               mexp,sign;
  56.  
  57. #ifdef national /* ordering of words in a floating points number */
  58.     int n0=1,n1=0;
  59. #else    /* national */
  60.     int n0=0,n1=1;
  61. #endif    /* national */
  62.  
  63.     mexp=px[n0]&0x7ff00000;
  64.     if(mexp==0x7ff00000) return(x); /* cbrt(NaN,INF) is itself */
  65.     if(x==0.0) return(x);        /* cbrt(0) is itself */
  66.  
  67.     sign=px[n0]&0x80000000; /* sign= sign(x) */
  68.     px[n0] ^= sign;        /* x=|x| */
  69.  
  70.  
  71.     /* rough cbrt to 5 bits */
  72.     if(mexp==0)         /* subnormal number */
  73.       {pt[n0]=0x43500000;     /* set t= 2**54 */
  74.        t*=x; pt[n0]=pt[n0]/3+B2;
  75.       }
  76.     else
  77.       pt[n0]=px[n0]/3+B1;    
  78.  
  79.  
  80.     /* new cbrt to 23 bits, may be implemented in single precision */
  81.     r=t*t/x;
  82.     s=C+r*t;
  83.     t*=G+F/(s+E+D/s);    
  84.  
  85.     /* chopped to 20 bits and make it larger than cbrt(x) */ 
  86.     pt[n1]=0; pt[n0]+=0x00000001;
  87.  
  88.  
  89.     /* one step newton iteration to 53 bits with error less than 0.667 ulps */
  90.     s=t*t;        /* t*t is exact */
  91.     r=x/s;
  92.     w=t+t;
  93.     r=(r-t)/(w+r);    /* r-t is exact */
  94.     t=t+t*r;
  95.  
  96.  
  97.     /* retore the sign bit */
  98.     pt[n0] |= sign;
  99.     return(t);
  100. }
  101. #endif
  102.